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Abstract 

We study the production of charmonia and charmed hadrons in nucleus-nucleus 
collisions at SPS, RHIC, and LHC energies within the framework of the statistical 
hadronization model. Results from this model are compared to the measured yields 
and centrality dependence of J/ip production at SPS and RHIC energies. We explore 
and contrast the centrality dependence of the production of mesons with open and 
hidden charm at SPS, RHIC and LHC. The sensitivity of the results to various 
input parameters is analyzed in detail for RHIC energy. 
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1 Introduction 



The concept of a possible thermal origin of charmonia was initially introduced [1] to 
explain the measurement of the J/^/hadron ratio in nuclear collisions. Because of their 
large mass and small production cross section at thermal energies charm quarks are, 
however, not likely to be thermally produced. On the other hand, significant production 
of charm quark pairs takes place in initial hard collisions. This led to the idea of statistical 
hadronization of charm quarks [2,3], which has sparked an intense activity in this field 
[4,5,6,7,8]. At about the same time an independent effort [9,10] based on a kinetic model 
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has been developed to study charm production in heavy ion collisions at collider energies. 
Initial interest focussed on the available SPS data on J/ip production in Pb-Pb collisions. 
As we stress below, these data can indeed be described within the statistical approach, 
but only when assuming the charm production cross section to be enhanced beyond the 
perturbative QCD (pQCD) predictions. 

The largest differences between the results of the statistical coalescence scenario and 
earlier models are expected at collider energies. For example, for central Au-Au collision 
at RHIC energy, the Satz-Matsui approach [11] predicts a very strong (up to a factor 
of 20 [12]) suppression of J/ip yields as compared to a direct production. In the present 
approach this suppression is overcome by statistical recombination of J/ip mesons from 
the initially produced cc pairs, so that much larger yields are expected. 1 We therefore 
focus in this letter on the predictions of the model for open and hidden charm meson 
yields at RHIC (Au-Au) and LHC (Pb-Pb) energies, with particular emphasis on the 
centrality dependence of rapidity densities at RHIC. We, furthermore, present predictions 
for the energy dependence of the production of hadrons with open and hidden charm. 



2 Model description and input parameters 

The model assumes that all charm quarks are produced in primary hard collisions and 
equilibrate 2 in the quark-gluon plasma (QGP). In particular, no J/ip mesons are pre- 
formed in the QGP, implying complete color screening [11]. There are indications from 
recent lattice QCD calculations that, already at SPS energies, J/ip mesons can be disso- 
ciated in a deconfined medium via collisions with thermal gluons [13,14]. As noted earlier 
[2], the cross sections of charmed hadrons are too small to allow for their chemical equi- 
libration in a hadronic gas. We assume that open and hidden charm hadrons are formed 
at chemical freeze-out according to statistical laws. Consistent with the fact that, at SPS 
energy and beyond, chemical freeze-out appears to be at the phase boundary [15], our 
model implies that a QGP phase was a stage in the evolution of the fireball. The analysis 
of J/ip spectra at SPS [16] lends further support to the statistical hadronization picture 
where J/ip decouples at chemical freeze-out. A recent analysis of electron spectra at RHIC 
[17] also strenghtens the case for an early thermalization of heavy quarks. However, in 
that analysis it was pointed out that both the hydrodynamical approach and PYTHIA 
reproduce the measured single-electron spectra. As suggested meanwhile [18], the mea- 
surement of the elliptic flow (t> 2 ) of open charm hadrons and J/ip may help to disentangle 
the two scenarios. Recent work [19] in extracting the quark t> 2 from the one of hadrons 
may help towards this goal. Most important though, in our view, is a high-precision di- 
rect measurement of open charm. As expected, and as seen in ref. [17], the p t spectra 
predicted by PYTHIA and hydro are different in detail at low p t and differ manifestly at 
high p t (p t 3> mass of charm quark). On the other hand, a recent theoretical analysis [20] 

1 The absence of J/ip suppression at RHIC energy is also predicted by the kinetic model [9,10]. 

2 This implies thermal, but not chemical equilibrium for charm quarks. 
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indicates that charm quarks might not thermalize quickly because of their large mass. 



In statistical models charm production needs generally to be treated within the framework 
of canonical thermodynamics [21]. Thus, the charm balance equation required [2,3] during 
hadronization is expressed as: 

N dir = -a N^ hjgcKc) , 2 N th f i) 



Here N^ r is the number of directly produced cc pairs and /„ are modified Bessel functions. 
In the fireball of volume V the total number of open = n^V and hidden = n%V 
charm hadrons are computed from their grand-canonical densities n*^ and n'g, respectively. 
The densities of different particle species in the grand canonical ensemble are calculated 
following the statistical model introduced in [23]. All known charmed mesons and hyperons 
and their decays [22] have to be included in the calculations. 

The balance equation (1) defines a fugacity parameter g c that accounts for deviations of 
charm multiplicity from the value that is expected in complete chemical equilibrium. The 
yield of open charm mesons and hyperons % and of charmonia j is obtained from: 

w< = ^iS§ a,ld N ' = 9 W- (2) 



The above model for charm production and hadronization can be only used [2,27] if the 
number of participating nucleons N part is sufficiently large. Taking into account the mea- 
sured dependence of the relative yield of if)' to J/ip on centrality [25] in Pb-Pb collisions 
at SPS energy the model appears appropriate [2,27] if N part >100. 

To calculate the yields of open and hidden charm hadrons for a given centrality and 
collision energy one needs to fix a set of parameters in Eq.(l) and (2): 

i) A constant temperature of 170 MeV and a baryonic chemical potential according 
to the parametrization ji b = 1270 [MeV] /(l + y/s^N [GeV]/4.3) [37] are used for our cal- 
culations. For SPS and RHIC energies these thermal parameters are consistent, within 
statistical errors, with those required to describe experimental data on different hadron 
yields [23,24]. In Table 1 we present the values of these input parameters and the result- 
ing chemical potentials for isospin, strangeness and charmness for the canonical ensemble 
calculations at SPS, RHIC and LHC. 

ii) The volume of the fireball corresponding to a slice of one rapidity unit at midrapidity, 
Va v =i, is obtained from the charged particle rapidity density dN ch /dy, via the relation 
dN c h/dy = ji^Vaj/=i, where n'^ is the charged particle density computed within the 
thermal model. The charged particle rapidity densities (and total yields in case of SPS) 
are taken from experiments at SPS [28] and RHIC [29] and extrapolated to LHC energy. 
In Fig. 1 we show a compilation of experimental data on dN c h/dy at midrapidity for 
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Table 1 

Input values of T and and output chemical potentials for isospin, strangeness and charmness 



\J snn (GeV) 


17.3 


200 


5500 


T (MeV) 


170 


170 


170 


Hb (MeV) 


253 


27 


1 


M/3 (MeV) 


-10.1 


-0.96 


-0.037 


M5 (MeV) 


-68.6 


6.70 


0.304 


(MeV) 


-41.8 


-3.92 


-0.156 



central collisions in the energy range from AGS [30] up to RHIC 3 . The continuous line in 
Fig. 1 is the prediction [33] of the saturation model, the dashed line is a natter power law 
dependence, both arbitrarily normalized. We use the (^snn) ' 3 dependence of dN c h/dy to 
get the value of dN ch /dy ~2000 for LHC energy. Throughout this paper we define central 
collisions by the average value N part =350, which roughly corresponds to the topmost 5% 
of the nuclear cross section. For the centrality dependences we assume that the volume of 
the fireball is proportional to N part . 

iii) The yield of open charm N^ r is taken from next-to-leading order (NLO) pQCD cal- 
culations for p-p collisions [34] and scaled to nucleus-nucleus collision via the nuclear 
overlap function [35,36] Taa- For a given centrality: N^ r (N part ) = <j(pp — > cc)TAA{N par t)- 
The pQCD calculations with the MRST HO parton distribution function (PDF) are used 
here [34]. We note that the pQCD results, in particular at LHC energy, are sensitive on 
the choice of PDF and/or charm quark mass value [34]. 

The input values dN c h/dy and dN^/dy and the corresponding volume at midrapidity 
and enhancement factor are summarized in Table 2 for model calculations for different 
collision energies. 



Table 2 

Input (dN^/dy and dN^/dy) and output (VA y =i and g c ) parameters for model calculations 
at top SPS, RHIC and LHC for central collisions. 



yJs NN (GeV) 


17.3 


200 


5500 


dN ch /dy 


430 


730 


2000 


dN^/dy 


0.064 


1.92 


16.8 


V Ay=1 (fm 3 ) 


861 


1663 


4564 


9c 


1.86 


8.33 


23.2 



A constant Jacobian of 1.1 has been used to convert the dA^/dr/ data to dN^/dy. 
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Fig. 1. Experimental values of charged particle yields at midrapidity as a function of collisions 
energy for central collisions. The data points are from experiments at AGS (E877 [30]), SPS 
(NA49 [28], NA50 [31]) and RHIC (PHOBOS [29], BRAHMS [32]). The continuous and dashed 
lines are power law dependences (see text), the dotted line marks the full LHC energy for Pb-Pb 
collisions (y/sNN=5.5 TeV). 



3 Model results and predictions 



We first compare predictions of the model to 47r-integrated J/ip data at the SPS energy. 
The NA50 data [38,39] are replotted in Fig. 2 as outlined in [3]. The model results shown 
in Fig. 2 are going beyond that of [3] since we have included a complete set of charmed 
mesons and baryons as well as updated the values of the charm production cross section 
and the volume of the fireball. Recent results [28] of the NA49 Collaboration on total 
charged particle multiplicity in central collisions yield N c h = 1533, leading to a fireball 
total volume V=3070 fm 3 and to a corresponding total yield of thermal open charm pairs 
of A^=0.98. This is to be contrasted with A^"" =0.137 from NLO calculations [34], leading 
to a value of g c =0.78. Although g c is here close to unity, this obviously does not indicate 
that charm production appears at chemical equilibrium, as the suppression factor is a 
strongly varying function of the collision energy. We have already indicated that, within 
the time scales available in heavy ion collisions, the chemical equilibration of charm is 
very unlikely both in confined and deconfined media [40]. 

In Fig. 2 we show the comparison between the results of our model and NA50 data for two 
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Fig. 2. The centrality dependence of J/ip production at SPS. Model predictions are compared 
to 47r-integrated NA50 data [38,39]. Two curves for the model correspond to the values of 
iV* r from NLO calculations (dashed line) and scaled up by a factor of 2.8 (continuous line). 
The dash-dotted curve is obtained when considering the possible NA50 A^ ar t-dependent charm 
enhancement over their extracted p— p cross section [41] (see text). 



different values of N^ r : from NLO calculations [34] and scaled up by a factor of 2.8. Using 
the NLO cross sections for charm production scaled by the nuclear overlap function, the 
model understimates the measured yield. To explain the overall magnitude of the data, 
one needs to increase the N^ r yield by a factor of 2.8 as compared to NLO calculations. 
We mention in this context that the observed [41] enhancement of the di-muon yield 
at intermediate masses has been interpreted as a possible indication for an anomalous 
increase of the charm production cross section [42]. A third calculation (resulting in the 
dash-dotted line in Fig. 2) is using the NLO cross section scaled-up by 1.6, which is the 
ratio of the open charm cross section estimated by NA50 for p-p collisions at 450 GeV/c 
[41] and the present NLO values. For this case the N part scaling is not the overlap function, 
but is taken according to the measured di-muon enhancement as a function of N part [41] . 
The resulting J / ip yields from the statistical model are on average in agreement with the 
data, albeit with a natter centrality dependence than by using the nuclear overlap function. 
Thus our charm enhancement factor of 2.8 needed to explain the J/ip data is very similar 
to the factor needed to explain the intermediate mass dilepton enhancement assuming 
that it arises exclusively from charm enhancement [41]. We note, however, that other 
plausible explanations exist of the observed enhancement in terms of thermal radiation 
[43,44]. 
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The drop of the J/ip yield per participant observed in the data for N part >350 (see 
Fig. 2) is currently understood in terms of further J /if) dissociation due to energy density 
fluctuations for a given overlap geometry [45] . Alternatively, the drop has been connected 
with a possible trigger bias in the data [46]. Note that this second drop is only visible 
in the minimum bias analysis [38]. These aspects are not included in the current model, 
however when being incorporated the J /if) yield for very central collisions can be as well 
reproduced [7,8]. 

Table 3 

Mid-rapidity densities for open and hidden charm hadrons, calculated for central collisions at 
SPS, RHIC and LHC. 



\/ s NN l" ev ) 


1 7 3 




oow 


dN D+ /dy 


0.010 


0.404 


3.56 


dN D -/dy 


0.016 


0.420 


3.53 


dN D o/dy 


0.022 


0.888 


7.80 


dN^o/dy 


0.035 


0.928 


7.82 


dN Ds+ /dy 


0.012 


0.349 


2.96 


dN Ds -/dy 


0.009 


0.338 


2.95 


dN A Jdy 


0.014 


0.153 


1.16 




0.0012 


0.117 


1.15 


dNj/^/dy 


2.55-10~ 4 


0.011 


0.226 


dAV/dy 


0.95-10- 5 


3.9740" 4 


8.46-10- 3 



We turn now to discuss our model predictions for charmonia and open charm production 
at collider energies and compare them with the results obtained at SPS. Notice that from 
now on we focus on rapidity densities, which are the relevant observables at the colliders. 
In Table 3 we summarize the yields for a selection of hadrons with open and hidden charm. 
All predicted yields increase strongly with beam energy, reflecting the increased charm 
cross section and the concomitant importance of statistical recombination. Also, ratios of 
open charm hadrons evolve with increasing energy, reflecting the corresponding decrease 
in the charm chemical potential. 

Model predictions for the centrality dependence of J/ip and D + rapidity densities nor- 
malized to Np ar t are shown in Fig. 3. The results for J /if) mesons exhibit, in addition 
to the dramatic change in magnitude, a striking change in the shape of the centrality 
dependence. In terms of the model this change is a consequence of the transition from the 
canonical to the grand-canonical regime [3]. For D + -mesons, the expected approximate 
scaling of the ratio D + /N part oc N^t (dashed lines in Fig. 3) is only roughly fulfilled due 
to departures of the nuclear overlap function from the simple N pa l t dependence. 

The results summarized in Table 3 and shown in Fig. 3 obviously depend on two input 
parameters, dN c h/dy and dN^ r /dy. For LHC energy, neither one of these parameters is 
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Fig. 3. Centrality dependence of rapidity densities of J /if) (upper panel) and D + (lower panel) 
mesons per participant at SPS, RHIC and LHC. Note the scale factors for RHIC and SPS ener- 

1/3 

gies. The dashed lines in the lower panel represent N art dependences normalized for N par t=350. 

well known. An increase of charged particle multiplicities by up to a factor of three beyond 
our "nominal" value diV c / l /d?/=2000 for central collisions is conceivable. The saturation 
model, e.g., leads to a prediction of [33] dN c h/dy ~2300. However, due to quite large 
uncertainties on the amount of shadowing at LHC energy, these results may be still 
modified. The yield of dN^ r /dy is also not well known at LHC energy. Although these 
uncertainties affect considerably the magnitude of the predicted yields, their centrality 
dependence remains qualitatively unchanged: the yields per participant are increasing 
functions of N part . We also note here that, while detailed predictions differ significantly, 
qualitatively similar results (see ref. [47]) have been obtained for a kinetic model study of 
J /if) production at the LHC. 



In Fig. 4 we present the predicted centrality dependence of the J/ ip rapidity density nor- 
malized to N part for RHIC energy ( v /J/v r /v=200 GeV). The three panels show its sensitivity 
on dN ch /dy, dN^ r /dy, and (freeze-out) temperature T. The calculations are compared 
to experimental results of the PHENIX Collaboration [48]. Note that the point for the 
central collisions is the upper limit extracted by PHENIX for 90% C.L. [48]. The experi- 
mental data have been rescaled according to our procedure [36] to calculate N part and the 
number of binary collisions, N co u. Using the centrality intervals of PHENIX, we obtain 
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Fig. 4. Centrality dependence of rapidity densities of J/ifi mesons at RHIC. Upper panel: sen- 
sitivity to dN c h/dy; middle panel: sensitivity to dN^ r /dy; lower panel: sensitivity to T. The 
calculations are represented by lines. The dots are experimental data from the PHENIX col- 
laboration [48]. Note that the point for the central collisions is the upper limit extracted by 
PHENIX for 90% C.L. [48]. 

N co u = 275 and 782 (the differences to the PHENIX values of 296 and 779, respectively, 
are small). We use (X7VAr=42 mb and a Woods-Saxon density profile. 

Within the still large experimental error bars, the measurements agree with our model 
predictions. In Fig. 4 only the statistical errors of the mid-central data point are plotted. 
The systematic errors are also large [48]. peripheral collisions. A stringent test of the 
present model can only be made when high statistics J/ip data are available. In any case, 
very large J /ip suppression factors as predicted, e.g., by [12] seem to be not supported by 
the data. 
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We turn now to a more detailed discussion of the sensitivity of our calculations to the 
various input parameters as quantified in Fig. 4. First we consider the influence of a 
10% variation of dN ch /dy on the centrality dependence of J/ip yield. Note that the total 
experimental uncertainty of dN ch /dr] (which is for the moment the measured observable 
for most experiments) at RHIC is below 10% [29,32]. The sensitivity on the dN ch /dy 
values stems from the volume into which the (fixed) initial number of charm quarks is 
distributed. The smaller the particle multiplicities and thus also the fireball volume, the 
more probable it is for charm quarks and antiquarks to combine and form quarkonia. 
That is why one sees, in the top panel in Fig. 4, that the J/ip yield is increasing with 
decreasing charge particle multiplicity. 

The sensitivity of the predicted J/ip yields on diV* r /dy is also straightforward. The larger 
this number is in a fixed volume the larger is the yield of charmed hadrons. In case of 
charmonia the dependence on dN^ r /dy is non-linear due to their double charm quark 
content, as reflected by the factor g 2 c in equation (2). To illustrate the sensitivity of the 
model predictions on dN^/dy, we exhibit the results of a 20% variation with respect to 
the value given in Table 2. The open charm cross section is not yet measured at RHIC. 
However, some indirect measurements can be well reproduced, within the experimental 
errors, by PYTHIA calculations using a p-p charm cross sections scaled with the number 
of collisions N co u of 650 /xb [50]. The corresponding value at v /JJv r /v=130 GeV is 330 /xb [49]. 
For comparison, the NLO pQCD values we are using [34] are 390 and 235 /xb, respectively. 
Despite the still large experimental uncertainties, this discrepancy, recognized earlier in 
ref. [10], needs to be understood. We note that, dependent on the input parameters used 
in the NLO calculations [34,10], possible variations of the open charm production cross 
section for the RHIC energy are of the order of ±20%. In terms of our model this variation 
corresponds to about a ±30% change in J / ip yield which is also centrality dependent (see 
middle panel in Fig. 4). If we use the PHENIX p-p cross section of 650 /xb, the calculated 
yield is a factor 2.5 larger for N part —350 and increases somewhat stronger with centrality, 
in disagreement with the data. As apparent in Fig. 4, the predictive power of our model 
or of any similar model relies heavily on the accurate knowledge of the overall charm 
production cross section. A simultaneous description of the centrality dependence of open 
charm together with J/ip production is, in this respect, mandatory to test the concept 
of the statistical origin of open and hidden charm hadrons in heavy ion collisions at 
relativistic energies. 

The apparent weak dependence of J/ip yield on freeze-out temperature, seen in Fig. 4, 
may be surprising. This is particularly the case since in the thermal model considered in [1] 
this dependence was shown to be very strong. In our model this result is a consequence of 
the charm balance equation (2). The temperature variation leads, obviously, to a different 
number of thermally produced charmed hadrons as in [1], but this is compensated by the 
g c factor. The approximate temperature dependence of g c and the J / ip yield are: 

2 c (T)~l/i\^~e^, %(T) = ^~e^ i (3) 
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Table 4 

Temperature dependence at RHIC, for central collisions. 



T 


160 


170 


180 


180/160 


dN% + /dy 


0.0368 


0.0554 


0.0775 


2.106 




0.672-10 4 


0.153-10 d 


0.311-10 6 


4.628 


9c 


12.85 


8.33 


5.77 


1/2.227 


&N D+ /dy 


0.415 


0.404 


0.391 


0.942 




0.0111 


0.0106 


0.0103 


0.928 




0.031 


0.037 


0.045 


1.452 



As a result of the small mass difference in the exponent the J/ ip yield exhibits only a weak 
sensitivity on T. This is in contrast to the purely thermal case where the yield scales with 
exp(— mj/^/T). The results are summarized in Table 4, where we present the number of 
J/t/j and D + mesons produced thermally {Nfj^ and N^ + , respectively) and within the 
frame of our statistical hadronization model. Also shown are the values of g c and the yield 
ratio ip'/J/ip for three different values of T. The compensations are evident in the last 
column of Table 4, where ratios for the two extreme temperatures are calculated. The only 
exception is the ratio ip'/J/ip, which is obviously identical in the statistical hadronization 
scenario and in the thermal model and coincides, for T~170 MeV, with the measured 
value at SPS [2]. 
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Fig. 5. Centrality dependence of rapidity densities of J/tp mesons at RHIC for different ra- 
pidity window sizes. The lines are calculations, the dots are experimental data from PHENIX 
collaboration [48] (the point for the central collisions is the upper limit for 90% C.L.). 
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Most of our results presented above are obtained considering a one unit rapidity window 
at midrapidity, while results for the full volume were presented only for the SPS. Unlike 
the kinetic model of Thews et al. [9,10], our model does not contain dynamical aspects 
of the coalescence process. However, in our approach, the width of the rapidity window 
does influence the results in the canonical regime. For the grand-canonical case, attained 
only at LHC energy, there is no dependence on the width of the rapidity window, due to 
a simple cancellation between the variation of the volume, proportional to the rapidity 
slice in case of a flat rapidity distribution, and the variation of N^ r , also proportional 
to the width of the rapidity slice. In Fig. 5 we present the centrality dependence of J/ip 
rapidity densities for RHIC energy and for different rapidity windows Ay from 0.5 to 3. 
The dependence on Ay resembles that of the kinetic model [10], but is less pronounced. 
The available data are not yet precise enough to rule out any of the scenarios considered. 
However, for the kinetic model [10], the cases of small Ay seem to be ruled out by the 
present PHENIX data. 

We stress in this context that the size of Ay window has a potentially large impact on 
the results at SPS energy. It is conceivable that no charm enhancement is needed to 
explain the data if one considers a sufficiently narrow rapidity window for the statistical 
hadronization. 

The results presented above were obtained under the assumption of statistical hadroni- 
zation of quarks and gluons. In addition, for charm quarks the yields were constrained 
by the charm balance equation that was formulated in the canonical ensemble. We have 
further assumed that charm quarks are entirely produced via primary hard scattering and 
thermalized in the QGP. No secondary production of charm in the initial and final state 
was included in our calculations. Final state effects like nuclear absorption of J/ip [51] are 
also neglected. 



4 Conclusions 



We have demonstrated that the statistical coalescence approach yields a good description 
of the measured centrality dependence of J/ip production at SPS energy, albeit with a 
charm cross section increased by a factor of 2.8 compared to current NLO perturbative 
QCD calculations. Rapidity densities for open and hidden charm mesons are predicted to 
increase strongly with energy, with striking changes in centrality dependence. First RHIC 
data on J / ip production support the current predictions, although the experimental errors 
are for the moment too large to allow firm conclusions. We emphasize that the predictive 
power of our model relies heavily on the accuracy of the charm cross section, not yet 
measured in heavy-ion collisions. The statistical coalescence implies travel of charm quarks 
over significant distances e.g. in a QGP. If the model predictions will describe consistently 
precision data this would be a clear signal for the presence of a deconfined phase. 
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